Transport properties controlled by a thermostat: An extended dissipative particle 
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o 
o 

(N 

Oh: 

m 



o 



I 

o 
o 



> 

p 

G\ 
O 
1^ 
O 



X 



Christoph Junghans, Matej PraprotnikQ and Kurt Kremer 
Max- Planck- Institut fiir Polymerforschung, Ackermannweg 10, D-55128 Mainz, Germany 

We introduce a variation of the dissipative particle dynamics (DPD) thermostat that aUows for 
controUing transport properties of molecular fluids. The standard DPD thermostat acts only on a 
relative velocity along the interatomic axis. Our extension includes the damping of the perpendicular 
components of the relative velocity, yet keeping the advantages of conserving Galilei invariance and 
within our error bar also hydrodynamics. This leads to a second friction parameter for tuning 
the transport properties of the system. Numerical simulations of a simple Lennard- Jones fluid and 
liquid water demonstrate a very sensitive behaviour of the transport properties, e.g., viscosity, on the 
strength of the new friction parameter. We envisage that the new thermostat will be very useful for 
the coarse-grained and adaptive resolution simulations of soft matter, where the diffusion constants 
and viscosity of the coarse-grained models are typically too high/low, respectively, compared to 
all-atom simulations. 
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I. INTRODUCTION 

Using the optimal set of degrees of freedom (DOFs) 
in computer simulations of soft matter guarantees effi- 
ciency, accuracy and avoids huge amounts of unnecessary 
detail, which might even obscure the underlying physics. 
This very idea is exploited in systematic coarse-graining 
efforts for modern computational materials science and 
biophysics problemsi^iS,, where atomistic simulations are 
usually beyond the possibilities of current and near fu- 
ture computers. A similar philosophy of reducing the 
number of DOFs is also employed by representing clus- 
ters of molecules with soft particles when simulating flu- 
ids on a mesoscopic scale using dissipative particle dy- 
namics (DPD)^i^i^ii^. However, in typical soft matter 
systems different time- and length-scales are intrinsically 
interconnected and a multiscale modeling approach is 
required to tackle such problems in the most efficient 
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Recently, we have proposed an adaptive resolution 
scheme (AdResS) that couples the atomistic and coarse- 
grained levels of detai l^''i-^^'-^^i^°'^-^i^^ . Due to the re- 
duction in DOFs upon coarse-graining, which elimi- 
nates the fluctuating forces associated with those missing 
molecular DOFs, the coarse-grained molecules typically 
move faster than the corresponding atomistically resolved 
ones^^i^. Although the accelerated dynamics is advan- 
tageous in some cases it can turn out to be problematic 
if one is really interested in dynamics in situations where 
two levels of resolutions are used within one simulation, 
as in the case of AdResS. To overcome this problem one 
can couple different classes of DOFs to the Langevin ther- 
mostat with different friction constants^i^iSliS^. We 
have shown in the example of liquid water that the 
coarse-grained dynamics can be slowed down by increas- 
ing the effective friction in the coarse-grained system^^. 
However, it is well known that the Langevin thermostat 
does not reproduce the correct hydrodynamics, i.e., the 
hydrodynamic interactions are unphysically screened. In 



order to correctly describe hydrodynamic interactions, 
one has to resort to the DPD thermostat^. 

In the past years DPD has established itself as a useful 
thermostat for soft matter simulations. The DPD ther- 
mostat is known to have several good properties, i.e., the 
thermostat satisfies Newton's third law by construction 
and owing to mass, momentum and temperature conser- 
vation, hydrodynamics is also correctly reproduced'^ ^. As 
it turns out, however, the DPD thermostat in its stan- 
dard form is not capable of controlling liquid properties 
such as viscosity and diffusion constant^. The aim of this 
work is to extend the DPD equations in such a way that 
these quantities can be tuned by changing the parameters 
of the thermostat. We consider the most general version 
of DPD^ and exploit the terms, which are not used in the 
standard DPD approach, i.e., the damping of the per- 
pendicular components of the relative particle velocities. 
This allows us to tune the viscosity of the coarse-grained 
liquid to match that of the all-atom counterpart while 
conserving the virtues of the standard DPD thermostat. 
This should be very useful for both coarse-grained and 
adaptive resolution simulations of soft matter. 

The article is organized as follows: In Sec. II the stan- 
dard and new DPD thermostats are presented. The simu- 
lation setup is described in Sec. HI. The results of molec- 
ular dynamics (MD) simulations of the Lennard- Jones 
fiuid and liquid water are reported in Sec. IV, followed 
by the conclusions in Sec. V. 



II. THERMOSTATS 

A. Standard DPD thermostat 

Newton's equations of motion are used in microcanon- 
ical NVE MD simulations and generate dynamics with 
constant energy. In order to run simulations in the canon- 
ical NVT ensemble the equations of motion have to be 
modified. In Langevin dynamics two additional forces are 
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introduced, a damping and a random force, whose ratio 
defines the temperature. The DPD equations of motion 
(used by the DPD thermostat) are then given hy^ 



and 



Pi 



(1) 



(2) 



where Fp denotes the conservative force on the ith par- 
ticle. The damping and random forces, can be split up 
in particle pair forces as 



FD ^ V F^ 
where the dissipative force reads as^ 



F, 



D 



and the random force is given by 



(3) 



(4) 



(5) 



(6) 



Flj = (Tllu;^(ry-)e,,f,, . 

In these equations the relative velocity Vij = — Vj be- 
tween the ith and jth particle was introduced, while 
denotes the unit vector of the interatomic 
fj. C" is the friction constant and (jH the noise strength. 
w^{rij) and w^(ry) are r-dependent weight functions. 
These are connected by the fluctuation-dissipation theo- 
rem (see Eqn. ©)■ The variable Qij is symmetric in the 
particle indices (O^ — Qji) and has the following first 



(e.,(<)) = o 



(7) 



and second moment 



{e,j{t)eki{t')) = 2{5,kSji + SaSjk)Sit - t') . (8) 
The fluctuation-dissipation theorem^ reads as 

(cTii)2^kBrc", (9) 



and 



(w"(r))' = w°(r). 



(10) 



The above DPD equations conserve the total momentum 
and reproduce correctly the hydrodynamics interactions 
in the system. However, previous studies^ii? have shown 
that the strength of the friction C'l does not influence the 
viscosity in linear order. In order to be able to tune the 
value of the viscosity of the system while preserving all 
of the virtues of the standard DPD thermostat presented 
briefly above we introduce in the next subsection its ex- 
tension named "Transverse DPD thermostat" . 



B. Transverse DPD thermostat 

We generalize Eqn. ([5]) and Eqn. ([6]) as 



and 



(11) 



(12) 



where C and cr are the friction constant and the noise 
strength of the generalized thermostat, respectively (see 
the text below). P ij{rij) is a projection operator 



p ^ pT ^ p2 



(13) 



which is symmetric in the particle indices ( P ij = P ji)- 
The scalar noise (see Eqn. ([8])) is replaced by a noise 
vector 9ij 

(4-(*) ® ^kiit')) = 2Y{S,kS,i ~ SuS,k)S{t ~ t') , (14) 

which is antisymmetric in the particle indices {9ij — 
—Oji) due to the symmetry of the projection operator 
and the antisymmetry of the pair force (Newton's third 
law). The corresponding Fokker-Planck operator >C is a 
sum of two parts: the deterministic part Cu 



E 



_d_ dH^_d_ m 

dfi dpi dpi dfi 



(15) 



and the generalized DPD part >Cdpd 
d 



P. 



^ OPt 



P 



dn on 

dpi dpj 

d_ _ 
dpi 



(16) 



The equilibrium condition £e~''^ = (£D+'CDPD)e~'^^ — 
then yields the dissipation-fluctuation-thcorcm in the 
same form as given by Eqn. ([9]) and Eqn. (fTO| . 

For the case where we choose the projector along the 
interatomic axis between particle i and j 



P ij {'''ij ) ~ ^ij ® ^ij ' 



(17) 



we retain the standard DPD thermostat. 

Alternatively, one can project on the plane perpendic- 
ular to the interatomic axis 



P ij i^ij ) 



(18) 



The space defined by the projector pS]) is orthogonal to 
the case of the standard DPD and introduces an exten- 
sion of the DPD thermostat, i.e., the Transverse DPD 
thermostat. Note that owing to this orthogonality the 
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new thermostat can be used in combination with the 
standard one. This enables us to adjust at the same 
time two friction constants C" and (-^ for the standard 
and Transverse DPD thermostats, respectively. Galilei 
invariance remains valid by construction while hydrody- 
namics is conserved within our error bar—. 

Our basic assumption is that in contrast to the stan- 
dard DPD the viscosity is very sensitive to the damp- 
ing perpendicular to the interatomic axis. This damping 
mimics the shear of those DOFs that were integrated 
out in the coarse-graining procedure. In a system with 
two particles the stochastic forces of the Transverse DPD 
thermostat act in the same direction as the shear forces. 
The mean force acting on a particle in the sheared sys- 
tem with more than two particles is hence a sum of two 
contributions: a force coming from the Transverse DPD 
thermostat and another one originating from the shear- 
ing of the probe. Therefore, the shear viscosity in a sim- 
ulation with the Transverse DPD thermostat is always 
higher than with the standard one. In the Green-Kubo 
picture this additional viscosity arises from the projected 
velocity-velocity autocorrelation function, which is de- 
rived by the Mori-Zwanzig formalism^i^^i^. The exact 
derivation is beyond the scope of this paper and will be 
presented elsewhere. Here, we shall demonstrate this by 
the results of our numerical experiments presented in the 
results section. 



III. SIMULATION SETUP 



B. Liquid Water 

All-atom water NVT simulations at ambient condi- 
tions are performed using the rigid TIP3P water model^. 
The electrostatics is described by the reaction field (RF) 
method, in which all molecules outside a spherical cav- 
ity of a molecular based cutoff radius i?c = 9 A are 
treated as a dielectric continuum with a dielectric con- 
stant crf = 80^1^1^. Typically, all-atom water sim- 
ulations are carried out using global thermostats, e.g., 
BerendseniS, Nose-Hoover'*'^'''^, Nose- Hoover chains^ 
thermostats, that dissipate the energy uniformly in the 
system. Local thermostats, e.g., Langevin, DPD^, 
Andersen^, Lowe-Andersen~2i^i^, Nose-Hoover-Lowe- 
Andersen^ thermostats, that dissipate energy on a spa- 
tially localized scale are usually used in coarse-grained 
simulations. Here, in order to reproduce the hydrody- 
namics correctly we employed the DPD thermostats^ 
with the friction constant — 0.5ps~^ and cutoff radius 
Re- The constant C" is small compared to the intrin- 
sic friction coefficient ^ of the TIP3P water system, i.e., 
^ = 288.6ps~^, so that the stochastic dynamics yields 
the correct dynamic o^^i^^ . For the coarse-grained wa- 
ter simulations with the Transverse DPD thermostat we 
have employed the single-site water model from Rcf.^^<, 
which reproduces essential thermodynamics and struc- 
tural properties, e.g., the pressure, density, and radial 
distribution functions, of the all-atom rigid TIP3P water 
model at standard conditions. Other simulation details 
are the same as given in Refi^. 



All simulations of the Lennard-Jones (LJ) liquid 
and liquid water are performed using the ESPResSo 
package'^i. 



IV. RESULTS 



Lennard-Jones Fluid 



A. Lennard-Jones Fluid 

We use the repulsive Weeks- Chandler- Andersen 
(WCA) potential 



(7Lj(r) = 4e 



^12 



a 

7^6 



(19) 



with the cutoff at Tc — 2^/®ct, a and e being the standard 
LJ parameters of length and energy. 

We chose as the weight function (see Eqn. (fTO|) ) for 
both thermostats a step function 



1, r <rc 
0, r > Tc 



(20) 



The simulations are carried out with a system consist- 
ing of 1000, 2000 and 4000 LJ particles at a temperature 
T = 1.2e/fcB and density p = Np^,t/V = l/(1.05a)3 in 
a cubic box with periodic boundary conditions. 



First, we checked in an equilibrium simulation the de- 
pendency of pressure and temperature on the strength of 
the friction. We set the reference temperature to 1.2e/kB 
and measured the instantaneous temperature defined as 



2i?kin 



3iV, 



part 



(21) 



where i?kin and iVpait are the kinetic energy and the num- 
ber of particles of the system, respectively. The relative 
deviation between the measured and reference tempera- 
ture was smaller then 1.2% for all strengths of friction 
and all combination of thermostats. The mean pressure 
at that temperature turned out to be 9.8±0.2e/a''^, which 
is in a perfect agreement with the results of previous 
studies^^. 

Next, we studied the dependency of the liquid trans- 
port properties, i.e., the diffusion constant and shear vis- 
cosity, on the friction constants C" and for the stan- 
dard and Transverse DPD thermostats, respectively. 
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FIG. 1; Diffusion constant (4000 LJ particles) as a function of 
the friction measured in equilibrium for different thermostats. 
In the case of the combined DPD thermostat only the strength 
of the friction parameter i^"'" was varied while the friction for 
the standard one was held constant at the value i^" = 1.0. 
The errors were obtained by averaging over several runs and 
Jackknife analysisS'^. 



1. Diffusion constant 

The diffusion constant was computed from the particle 
displacements using the Einstein relation 



D 



lim 

t — >oo 



\m-m\ 

6t 



(22) 



A small influence on this constant from the standard 
DPD thermostat'^*' is expected, but a considerable one 
from the new Transverse DPD thermostat. Former 
results^" for the standard DPD thermostat could be con- 
firmed. The value of the diffusion constant approaches 
the equilibrium value {D = 0.08a^/r) for vanishing fric- 
tion of the Transverse DPD thermostat (see Fig. [1]). The 
diffusion constant D is very sensitive on the friction (^-^ 
for the Transverse DPD thermostat. By changing (-^ it 
is therefore possible to tune the diffusion constant. 



2. Shear viscosity 

The viscosity was measured in nonequilibrium molecu- 
lar dynamics (NEMD) simulation by shearing the system 
with a constant shear rate in the y-direction2£ 



7 



dux 
dy 



(23) 



The viscosity can then be determined by the following 
simple formula 



77 = 



7^2 



(24) 



FIG. 2: Shear viscosity measured with the NEMD algorithm 
for different thermostats (4000 LJ particles) and the different 
shear rates (0.1 and 0.01 shown), which then are extrapolated 
to vanishing shear rate. The errors are obtained by Jackknife 
analysis^. 



where F is the mean force (momentum transfer per time 
unit). The apparent shear viscosity rj measured in NEMD 
simulation approaches the equilibrium viscosity with de- 
creasing shear rate. 

We found nearly no dependency on the strength of the 
friction C" for the standard DPD thermostat, (see also 
Ref.'^^). In contrast the friction (-^ for the Transverse 
DPD thermostat gives a very sensitive means of control- 
ling the viscosity (see Fig. [5]). In the case of the combi- 
nation of both thermostats the shear viscosity is mostly 
controlled by the Transverse DPD thermostat. In the 
limit of a vanishing shear rate a value of 2.45±0.07er/(T"^ 
was extrapolated, which matches former results^. This 
extrapolation has to be done due to short characteris- 
tic timescale of the LJ system. Additionally, we also 
checked that the equilibrium correlation of the pressure 
tensor is in accordance with the Grcen-Kubo picture: all 
(non-auto) off-diagonal - (off-)diagonal elements are un- 
correlated, (60 of 81 possible elements are zero). 

For higher values of the apparent viscosity becomes 
increasingly dependent on the shear rate. This is due to 
the fact that the dynamics is controlled in this regime by 
the thermostat forces (that are linear in C^) and hence 
we end up measuring the "viscosity of the thermostat" . 



B. Tuning the Dynamics of Water 

Having shown that the new thermostat enables us to 
tune the diffusion constant and viscosity of a simple liq- 
uid we apply it to an important physical example, i.e., 
liquid water at ambient conditions. 

We first check that the structural properties do not 
depend on the thermostat and also that we obtain con- 
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duce both the structure and the dynamics of the all-atom 
liquid water with the single-site coarse-grained water 
model. This is essential for synchronizing the timescales 
of the all-atom and coarse-grained regimes in the adap- 
tive resolution MD simulations^^. 
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sistency between the coarse-grained and atomistic simu- 
lations. The center-of-mass radial distribution functions 
of the all-atom and coarse-grained system using different 
values of match within the line thickness (see Fig. [3]). 

There is an intrinsic timescale difference in the diffusive 
dynamics of the coarse-grained water system because of 
the reduced number of DOFs, i.e., the self-diffusion con- 
stant for the coarse-grained water model is approximately 
2 times larger then the corresponding all-atom one using 
the Langevin thermostat with the same background fric- 
tion in both cases^S. We used — Q.8ps~^ for the 
Transverse DPD thermostat (C" = 0) to match the dif- 
fusion constant of the coarse-grained water model to the 
corresponding value D — 3.0 • 10~^m^/s obtained from 
the all-atom simulation with the standard DPD thermo- 
stat. With a friction strength of — 0.6ps~^, which 
is as desired very close to the above value for matching 
the diffusion constants, we were also able to match the 
viscosity to the desired value 77 — 0.5 ± 0.1 • lO^^Pa • s^^ 
for the TIP3P water model (from our atomistic simula- 
tion with C'l — 0.5ps~^). In this case the characteristic 
times for the atomistic model are much longer than the 
time scale of the shearing (I/7). Therefore even with a 
shear rate of 7 = 0.01 we are in the no shear limit and 
hence no extrapolation is required. The obtained diffu- 
sion constants and viscosities are in good agreement with 
the pubHshed datai^. 

Thus, employing the new thermostat one can repro- 



FIG. 4: The mean square displacements over time plot of the 
all-atom ((^" = 0.5ps~^ , = 0) and several coarse-grained 
simulations ((^" = and (^l' = 0.5ps~^, (^i' = 0.875ps~^ and 



V. CONCLUSIONS 

In this paper we introduced an extension of the DPD 
thermostat that allows for controlling the transport prop- 
erties of molecular liquids, e.g., water, while preserving 
the hydrodynamics. 

The presented Galilean invariant local thermostat can 
be used in coarse-grained simulations to tunc the diffu- 
sion constant and viscosity of the system to the desired 
values. This opens up the possibility of reproducing the 
atomistic dynamics by coarse-grained simulations, as it 
is required, for example, in recently introduced adaptive 
resolution simulations. 
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